# -*- coding: utf-8 -*-
"""
Created on Wed Apr  2 22:32:26 2025

@author: Administrator
"""
# -*- coding: utf-8 -*-
"""
Created on Tue Apr  1 10:55:34 2025

@author: Administrator
"""
import numpy as np
import pandas as pd


class AerodynamicCalculator:
    def __init__(self, U, a, rho=1.225):
        """
        初始化 AerodynamicCalculator 类的实例。

        参数:
        U (float): 自由流速度
        a (float): 圆柱体半径
        rho (float, 可选): 流体密度，默认为 1.225 kg/m³
        """
        self.U = U
        self.a = a
        self.rho = rho

    def compute_pressure(self, z_samples):
        """
        计算给定复平面上点的压力系数。

        参数:
        z_samples (numpy.ndarray): 复平面上的点

        返回:
        numpy.ndarray: 压力系数数组
        """
        # 计算速度复势的导数，得到速度分量
        V = self.U * (1 - (self.a**2) / (z_samples**2))
        # 计算速度大小
        V_mag = np.abs(V)
        # 根据伯努利方程计算压力系数
        Cp = 1 - (V_mag / self.U)**2
        return Cp

    def export_report(self, filename="pressure_report.csv"):
        """
        生成压力系数的误差分析报告，并保存为 CSV 文件。

        参数:
        filename (str, 可选): 报告文件名，默认为 "pressure_report.csv"
        """
        # 在 0 到 2π 范围内生成 36 个等间距的角度值
        theta = np.linspace(0, 2*np.pi, 36)
        # 将角度转换为复平面上的点
        z_samples = self.a * np.exp(1j * theta)
        # 计算数值压力系数
        Cp_numerical = self.compute_pressure(z_samples)
        # 计算理论压力系数
        Cp_theory = 1 - 4 * np.sin(theta)**2  

        # 计算绝对误差
        absolute_error = Cp_numerical - Cp_theory
        # 计算相对误差
        relative_error = absolute_error / Cp_theory

        # 计算驻点（theta = 0°, 180°）的理论和数值压力系数
        theta_stagnation = np.array([0, np.pi])
        z_stagnation = self.a * np.exp(1j * theta_stagnation)
        Cp_numerical_stagnation = self.compute_pressure(z_stagnation)
        Cp_theory_stagnation = 1 - 4 * np.sin(theta_stagnation)**2

        # 计算极值点（theta = 90°, 270°）的理论和数值压力系数
        theta_extremum = np.array([np.pi/2, 3*np.pi/2])
        z_extremum = self.a * np.exp(1j * theta_extremum)
        Cp_numerical_extremum = self.compute_pressure(z_extremum)
        Cp_theory_extremum = 1 - 4 * np.sin(theta_extremum)**2

        # 创建包含误差分析结果的 DataFrame
        df = pd.DataFrame({
            'theta_deg': np.degrees(theta),
            'Cp_numerical': Cp_numerical,
            'Cp_theory': Cp_theory,
            'Absolute_Error': absolute_error,
            'Relative_Error': relative_error
        })
        # 将 DataFrame 保存为 CSV 文件
        df.to_csv(filename, index=False, float_format="%.4f")
        print(f"误差分析表已保存至 {filename}")

        # 输出驻点的验证结果
        print("驻点（theta = 0°, 180°）验证结果：")
        for i, theta in enumerate(theta_stagnation):
            print(f"theta = {np.degrees(theta):.0f}°:")
            print(f"  数值解 Cp: {Cp_numerical_stagnation[i]:.4f}")
            print(f"  理论解 Cp: {Cp_theory_stagnation[i]:.4f}")
            print(f"  绝对误差: {Cp_numerical_stagnation[i] - Cp_theory_stagnation[i]:.4f}")
            if Cp_theory_stagnation[i] != 0:
                print(f"  相对误差: {(Cp_numerical_stagnation[i] - Cp_theory_stagnation[i]) / Cp_theory_stagnation[i]:.4f}")
            else:
                print("  相对误差: 理论解为0，无法计算相对误差")

        # 输出极值点的验证结果
        print("极值点（theta = 90°, 270°）验证结果：")
        for i, theta in enumerate(theta_extremum):
            print(f"theta = {np.degrees(theta):.0f}°:")
            print(f"  数值解 Cp: {Cp_numerical_extremum[i]:.4f}")
            print(f"  理论解 Cp: {Cp_theory_extremum[i]:.4f}")
            print(f"  绝对误差: {Cp_numerical_extremum[i] - Cp_theory_extremum[i]:.4f}")
            if Cp_theory_extremum[i] != 0:
                print(f"  相对误差: {(Cp_numerical_extremum[i] - Cp_theory_extremum[i]) / Cp_theory_extremum[i]:.4f}")
            else:
                print("  相对误差: 理论解为0，无法计算相对误差")

    def parameter_sensitivity_test(self, U_values=[1.0, 5.0, 10.0], a_values=[0.5, 1.0, 2.0]):
        """
        进行参数敏感性测试，分析不同自由流速度和圆柱体半径对压力系数的影响。

        参数:
        U_values (list, 可选): 自由流速度列表，默认为 [1.0, 5.0, 10.0]
        a_values (list, 可选): 圆柱体半径列表，默认为 [0.5, 1.0, 2.0]

        返回:
        pandas.DataFrame: 包含参数和对应最大、最小压力系数的 DataFrame
        """
        results = []
        # 遍历不同的自由流速度和圆柱体半径组合
        for U in U_values:
            for a in a_values:
                # 创建新的 AerodynamicCalculator 实例
                calc = AerodynamicCalculator(U=U, a=a)
                # 在 0 到 2π 范围内生成 36 个等间距的角度值
                theta = np.linspace(0, 2*np.pi, 36)
                # 将角度转换为复平面上的点
                z_samples = a * np.exp(1j * theta)
                # 计算压力系数
                Cp = calc.compute_pressure(z_samples)
                # 记录参数和对应的最大、最小压力系数
                results.append({
                    'U': U,
                    'a': a,
                    'Max_Cp': np.max(Cp),
                    'Min_Cp': np.min(Cp)
                })
        # 创建包含结果的 DataFrame
        df = pd.DataFrame(results)
        # 将 DataFrame 保存为 CSV 文件
        df.to_csv('sensitivity_analysis.csv', index=False)
        print("参数敏感性分析表已生成")
        return df

# 测试代码
if __name__ == "__main__":
    # 创建 AerodynamicCalculator 实例
    calc = AerodynamicCalculator(U=5.0, a=1.0)
    # 生成误差分析报告
    calc.export_report()
    # 执行参数敏感性测试
    sensitivity_df = calc.parameter_sensitivity_test()
    print(sensitivity_df)